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A new least squares algorithm is proposed and investigated for fast frequency and 
phase acquisition of sinusoids in the presence of noise. This algorithm is a special case 
of more general adaptive parameter-estimation techniques . The advantages of the 
algorithms are their conceptual simplicity , flexibility and applicability to general situa- 
tions. For example , the frequency to be acquired can be time varying and the noise can 
be non-gaussian, non-stationary and colored. 

As the proposed algorithm can be made recursive in the number of observations , it 
is not necessary to have a priori knowledge of the received signal-to-noise ratio or to 
specify the measurement time. This would be required for batch processing techniques , 
such as the fast Fourier transform (FFT). The proposed algorithm improves the frequency 
estimate on a recursive basis as more and more observations are obtained. When the 
algorithm is applied in real time, it has the extra advantage that the observations need 
not be stored. The algorithm also yields a real time confidence measure as to the accuracy 
of the estimator. 


I. Introduction 

The problem of estimating the parameters of a sinusoidal 
signal has received considerable attention in the literature, see 
for example Refs. 1-7 and their references. Such a problem 
arises in diverse engineering situations such as carrier tracking 
for communications systems and the measurement of Doppler 
in position location, navigation and radar systems. 

A variety of techniques have been proposed in the literature 
to solve such problems including, to mention a few, the appli- 
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cation of the fast Fourier transform (FFT) (as in Refs. 1, 2), 
one and two dimensional Kalman filters based on a linearized 
model (Ref. 5), a modified extended Kalman filter that 
results in a phase locked loop (Ref. 6), or a digital phase locked 
loop derived on the basis of linear stochastic optimization 
(Ref. 7). 

The fact that there are so many different techniques to 
solve the problem indicates the importance of the problem. 
This, however, also implies that there is no single technique 
superior to all others in all possible situations and/or with 
respect to different criteria such as computational complexity, 
statistical efficiency, etc. 
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In this article we propose the application of the least 
squares parameter estimation technique to the estimation of 
an unknown frequency. The least squares algorithm has been 
extensively studied in the literature in terms of convergence, 
computational requirements, etc. (Refs. 8, 9), and has found 
varied applications in a wide variety of communication and 
signal processing problems. This is due to the relative simplicity 
of the least squares algorithm and its attractive convergence 
rates. In Ref. 9 for example, it has been shown that the 
algorithm exhibits an initial factorial convergence rate fol- 
lowed by exponential convergence. Such a convergence is 
very desirable in almost all estimation situations including the 
one under consideration. 

When the least squares (LS) algorithm is implemented via 
the fast algorithm of Ref. 10, the computational requirements 
of the algorithm compare favorably to the FFT algorithm. 
The least squares algorithm offers, in addition to the above 
discussed rapid initial convergence, several other desirable 
features. First the least squares algorithm provides final 
estimates of frequency, whereas FFT estimation requires use 
of a secondary algorithm to interpolate between frequencies. 
Secondly, using an exponentially weighted least squares algo- 
rithm, it is possible to track a time varying frequency. We 
compare the least squares algorithm to the FFT since the latter 
is “close” to the optimum in terms of the statistical 
efficiency (Ref. 1). 


II. The Signal Model 

Consider the problem of estimating an unknown frequency 
w d from the measurements y k , z k below 

y k = A sm (w d t k + <p) + n. k 

( 1 ) 

z k = A cos ( w d t k + <P) + n qk , k = l ,2, . . . 

Here the sequence , z k } represents the samples of the 
in-phase and quadrature components of a received signal 
s(t) obtained by demodulating s(/) by a carrier reference 
signal r(t ) and its 90-deg phase shifted version respectively, 
i.e., 

s(t) = A sin(w 0 f + 0 o ) + rt(O 

r(t) = 2 sin (vv. t + 0 c ), 0 = = ~ w o 

with n ik and n qk denoting the samples of the quadrature com- 
ponents of white noise n{t). The algorithm can be easily 
extended to the case where n{t) is a colored noise. 

With a power series expansion for the sine and cosine func- 
tions, the measurement equations can be written in alternative 
forms as follows: 


In Section II we present the signal model followed by the 
least squares algorithm in Section III. Section IV analyzes 
the estimation error of the algorithm. In Section V a few simu- 
lation examples are presented. The last section of the article 
contains some concluding remarks. 


y k = A sin (w d t k ) cos <p + A cos (w rf r fc ) sirup + n jk 
z k = A cos (w d t k ) cos <p - A sin (w d t k ) sirup + n qk 
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In the above approximation the terms of the order (w d t k ) n /n\ 
and smaller order have been ignored (assuming here that 
w d t k < «)• With obvious definitions, the measurement equa- 
tion can be written in a form “linear in parameters.” 


Z k =e'x k + n k (3) 

In the above, the prime (') denotes transpose, Z k = [y k z k \ , 
n k = [ n ik n qk\ > x k denotes the observable state vector [1 t k 
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t\ . . . /jj?' 1 ] and O' is the unknown parameter matrix. A stan- 
dard least square algorithm can be applied to estimate the 
unknown parameter matrix O' from the sequence of noisy 
observations Z k ,k =1,2 , . . . ,A r . 


III. Parameter Estimation via Least Squares 

The parameter matrix O' can be estimated by either a 
recursive or nonrecursive form. We consider in this article the 
nonrecursive form. The estimate of 6 on the basis of measure- 
ment Z k ,k= 1,2, . . . , N, denoted 0 N , is given by 



Each of the matrices XjXj and P” 1 is a Hankel matrix. That is, 
all the elements of each cross-diagonal are the same. The struc- 
ture of a Hankel matrix is very similar to that of a Toeplitz 
matrix wherein the elements along the various subdiagonals 
are equal. The fast algorithm of Ref. 10 for the solution of 
Toeplitz system of equations can be slightly modified so as 
to become applicable to the present problem. Thus, the com- 
putations in Eq. (4) can be made in order n(\og 2 n) 2 compu- 
tations, resulting in considerable reduction in the require- 
ment for large values of n. 

If the matrix inverse is precomputed, then with the algorith- 
mic properties of Ref. 10, the solution for can be obtained 
in approximately 6«log 2 « operations. In the implementations 
above, it is sufficient to store only the first row and column of 
PorP-i. 


where 0 < X < 1 is the exponential data weighting factor. One 
may refer to Refs. 8 and 11, for example, for an equivalent 
recursive update of d N . From On. the estimates of A , w d and 
<p can be obtained. 

A. Computational Requirements 

The algorithm of Eq. (4) requires an inverse of a symmetric 
(n X n) matrix once, requiring order n 2 computations. It may 
appear that the computation of each jc.jc' term requires n 2 
computations. However, detailed examination shows that only 
2 n computations are required. Thus, the total computations 
are equal to 6nN + 0(n 2 ). In practice, the matrix inverse can 
be precomputed, thus reducing the data dependent computa- 
tions to only 2nN + n 2 / 2. 

B. Fast Implementation of Least Squares Algorithm 

The matrix 


N 


P- 1 = 


1 - E 


/=! 


x.x. X N j 


in Eq. (4) has a very special structure as can be easily seen by 
explicit computation of the term XjXj of the summand. Thus, 
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C. Baseband Sampling 

In the case of baseband sampling, only the measurements 
{ y k } are available and the parameter matrix O' is of dimension 
n X 1 . In such an implementation, however, there may result 
an ambiguity of 7r radians in the phase estimate if the sign of 
w d is also unknown. 


IV. Estimation Error Analysis 

Assuming that the model Eq. (3) is exact (the dimension n 
of the parameter matrix in Eq. (2) is sufficiently high), then 
the substitution of Eq. (3) in Eq. (4) yields, 

° N = (t, x i x 'i j jz x i (x 'j 6 + n ? | 

( 5 ) 

A simple manipulation of Eq. (5) yields the estimation error 
Opf ~~ 0 — 0^ as 



As the^ state vector Xj is deterministic, and n } - is a zero mean 
process, 0 N has its mean equal to zero. The error covariance 
matrix of 0 N can also be evaluated in a straightforward manner. 
Post multiplying Eq. (6) by the transpose of and taking 
expected values of both sides. 
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x.x. X N " ; 


For the case of uniform sampling f ; . = jT s , where T s is the 
sampling period. Substituting for tj and letting T - N T s denote 
the observation period, 


X X x'itfw-n 
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i=i 


X 
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x i x i x 



( 7 ) 


Considering the case of X = 1 and recalling that fy} is a white 
noise sequence, 



(8) 


p, r ~2 ] _ ^ pr 6 1 1 

‘ vJ " 2 7V(jV + 1X2 N+ 1) T 7 A 7 

1 S 

In terms of the unsampled system, if the additive noise process 
has one-sided noise spectral density N 0 , then a 2 = 2NJT S . 
Thus, 

E l”Jl r ■ W 

where P = A 2 /2 is the received signal power and K has value 
approximately equal to 4 for low values of n. This is the same 
mean square error as for Maximum Likelihood estimation 
(Ref. 1; Ref. 12, Eq. 8.1 16). 


E [llriyll 2 ] = E [n 2 ] +E [n 2 . ] = a 2 

A. Frequency Estimation Error 

A simple approximate expression can also be obtained for 
the frequency estimation error when the amplitude A is known 
and uniform sampling is used. The frequency estimate w d can 
be obtained as 



We note here that in the derivation of Eq. (9), the approxi- 
mation error in Eq. (3) has been ignored. It is difficult to esti- 
mate the error due to such finite approximation. However, 
from a few computer simulations, it appears that for n > 
w d T = (w d T S )N, such error is small. 

B. Examples of Application to the DSN Receiver 

To keep the dimension n of the parameter matrix small, 
the following estimation method is proposed. Dividing both 
sides of Eq. (9) by wj max , and substituting T- n/w d max , one 
obtains 


When the amplitude A is also unknown, it can be replaced by 
its estimate given by, 


A N 



In the above expression, denotes the (/,/)th element of 
the parameter matrix 6 N . The error variance of these elements 
of interest is given by 

e [(^‘) ! ] t£ [(‘”) 2 ] aA: " 2 (2 »/) 1 

where K approaches a constant with the increase in the num- 
bers of observations N. 


For relatively small errors, the frequency estimation error 
n = w d ~ n ^ as variance approximately 


£[g}l 

W d ,max 





,max 


n 


3 


Selecting a value of 1/36 for the left hand side of the above 
equations allows one to express the maximum frequency 
uncertainty that can be resolved by the algorithm as a function 
of n. Thus 

« 3 p 

W d , max 216 N 

O 

The rationale for selecting the value of 1/36 for £[vvj]/ 
w d max is as Since the additive noise has Gaussian dis- 

tribution, one may assume that the frequency estimation error 
has Gaussian distribution with its standard deviation denoted 
by o ^ d . The above selection thus ensures that 3 o^ d < 

W d,maxl^- 



Example 1. For reception of Voyager 2 signals at DSS 13, a 
typical carrier power-to-noise spectral density ratio is 24.4 
dB-Hz. Let n = 8, and w d max = 652 rad/s. After an initial 
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estimation period of T = n w^j max , the receiver NCO fre- 
quency is adjusted by \v d . Thus, with an initial adjustment 
after T = 12.2 ms, the frequency offset is reduced to w d with 
o = 108.6 rad/s (17.4 Hz). Application of the algorithm for a 
subsequent period of 24.4 ms reduces the standard deviation 
to 6 Hz. In this manner, four applications of the algorithm 
bring down the standard deviation of the frequency offset to 
less than 0.7 Hz in a total time of 183 ms, from an initial 
frequency offset of 104 Hz. 

Example 2. If the initial uncertainty is only 20 Hz, then 
with a lower value of n equal to 5, after an estimation period 
of T = 40 ms, o~ d =18.4 rad/s (2.94 Hz). A frequency correc- 
tion at the end of this period and an estimation of the residual 
frequency offset for a period of 160 ms reduces of d to 0.36 Hz 
(. f d - w d /2ir). Thus with n = 5, the frequency uncertainty is 
reduced to Of d - 0.36 Hz in a total estimation period of 0.2 s. 

In an alternative approach to keep the value of n fixed and 
small with an increase in the total observation period, instead 
of resetting the frequency reference (making a correction in 
the NCO frequency), the time reference is reset to zero. Sub- 
sequent observations in Eq. (1) are now with respect to a 
different phase reference, say 0~. The application of the least 
squares algorithm to this set of observations then provides an 
estimate for 0 denoted 0l At this stage the observations in the 
second T s interval are processed to have a phase reference 0 
and are then combined with the first set of observations. 
Equivalently, it is required to post multiply the second sum on 
the right hand side of Eq. (4), obtained for the second sub- 
interval of T s, by the following matrix 


cos (A0) 
^sin (A0) 


-sin (A 0) 
cos (A0)_ 


A0 = 0-0 


the unknown frequency. To keep the computational burden of 
the simulations to a minimum, the dimension n was restricted 
to a small value and the observation period was also restricted 
to a small value. 

For frequencies much higher than one, the least squares 
algorithm, Eq. (4), was slightly modified. Thus, as w d max 
denotes an upper bound on the magnitude of unknown fre- 
quency, we define a normalized parameter matrix 6 by 6 = 
; / = 1, . . . , n\j - 1, 2. Defining a corresponding 
state vector x k by 


** = 


1 0 

0 w . 

cf, max 
0 0 


w 


0 

0 

n - 1 
d , max 


the measurement Eq. (3) may be rewritten as 

Z k = e'x k + n k (10) 


and the least squares algorithm can now be applied to estimate 
6 . The estimates of the elements of 0 are then obtained as 


e iJ = e iJ 

d y max 


Such a transformation leaves the previous error analysis 
invariant. However, for Finite dimensional approximation 
considered here, this makes the algorithm numerically more 
robust. The simulations for w d = 10 and w d = 100 are pre- 
cisely the same as in Figs. 1 through 5 with appropriate 
changes in scaling and are not presented separately. 


and add the result to the corresponding sum for the first 
T s interval. The first sum on the right hand side of Eq. (4) is 
simply multilied by a factor of 2. This procedure is extended 
in an appropriate manner to subsequent intervals, so as to 
obtain a final estimate for w d , and 0 based on the complete 
set of observations. 

V. Simulations 

Figures 1 through 5 present the frequency estimates obtained 
by the least squares algorithm. To avoid singularity of the 
matrix P ~ x , it was modified by the addition of a diagonal 
matrix el with e = 0.001. For convenience, the unknown fre- 
quency w d is taken to be 1 rad/s. From the simulations it is 
apparent that the frequency estimate comes close to the true 
frequency in a time equal to a fraction of the time period of 


In Figs. 1 through 5, 0 and £ represent the first and second 
row of the parameter matrix respectively. Thus in baseband 
sampling, only the 0 vector is estimated while in quadrature 
sampling, the estimates of both 0 and £ parameter vectors are 
available. From the figures it is apparent that the dimension n 
of the state vector lx k in the model Eq. (3) is approximately 
equal to w d T where T is the sampling interval. Due to over 
parameterization involved in the problem , there is a considerable 
amount of flexibility in the estimation of A, w d , 0 from the 
estimates of the elements of 6. Thus whereas in Figs. 1 
through 3, the amplitude A is assumed known, Figs. 4 and 5 
involve unknown A. A different order of computation can 
provide an estimate of A when baseband sampling is used. 
Here, we have reported results only for the frequency esti- 
mates; the phase esimates also converge at a fast conver- 
gence rate. 
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VI. Comparison with FFT Techniques 

An alternative technique for the fast frequency acquisition 
is via fast Fourier transform of sampled data. We observe that 
for the case of infinite observation time, both procedures are 
optimum and thus are equivalent. However, for finite observa- 
tion period 7 , the FFT has the limitation that the frequency 
estimates are quantized to intervals of 1/7 Hz. In the finite 
dimensional approximation of the LS algorithm this is not the 
case and sufficiently accurate estimates can be obtained by 
choosing n sufficiently large (finer sampling) even for low 
values of 7. 

The price for such an improvement is increased computa- 
tional requirements which is of order n\og 2 n (though higher 
than for FFT) if the matrix P is precomputed and is of order 
n\ogn log n if P must be computed on line. With the applica- 
tion of fast algorithms, the storage requirement of P is only 
In (not n 2 j 2). 

Also, note that the computational requirements here are 
dominantly decided by w d T and not by the number of sam- 
ples as is the case with FFT. 

It may also be mentioned that with the FFT algorithm 
there also exists a finite probability of the occurrence of an 
outlier (Ref. 1) and this causes a component of the frequency 
estimation error with a uniform probability density function 
over the complete frequency range of the FFT algorithm. As 


against this, the frequency estimation error with LS algorithm 
has a Gaussian distribution. 


VII. Conclusion 

This article has presented a fast algorithm based on the 
least squares parameter estimation technique. In Ref. 9 it is 
shown that the least squares algorithm exhibits a convergence 
phase wherein the convergence rate is factorial (the estimation 
error goes to 0 as \/k\ where k is the number of observations) 
followed by an exponential convergence rate. Our simulations 
also exhibit the same rapid initial convergence rates. Here of 
course, the estimation error does not approach zero because of 
a finite and low dimensional truncation of the model. From 
another viewpoint the algorithm may be perceived as a time 
domain dual of the FFT algorithm. Whereas the FFT algo- 
rithm transforms the data into frequency domain for the 
estimation/detection purpose, here the estimation is done 
directly in the time domain. This latter approach has several 
advantages. First by choosing X < 1 , it is possible to track the 
time varying frequency by recursive update techniques 
(Ref. 8). Moreover unlike the case of the FFT algorithm, the 
frequency estimates are not quantized to intervals of 1/7 Hz, 
which would be large for small observation interval 7. The 
price for these desirable features is in terms of increased com- 
putational requirement which in fast implementation of the 
algorithm could be of order n\ogn or n\ogn log n (depending 
upon the specific implementations), where n is approximately 
equal to w d 7, the product of the frequency uncertainty and 
the observation period. 
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FREQUENCY ESTIMATE 



DIMENSION = 6 (ESTIMATE ONLY*) 
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Fig. 1. Least squares algorithm: noise free case and 
baseband sampling, n = 6 
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DIMENSION = 4 (ESTIMATE ONLY* ) 

w d =iV' a = i 

X - 0.99 
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Fig. 3. Least squares algorithm: noise free case and 
baseband sampling, n = 4 


DIMENSION = 5 (ESTIMATE ONLY*) 


|9 2 /|1 - (t,) 2 l ,/2 l 


d - 2 1 v T i ' ' i 

X = 1 

SAMPLING RATE ~ 25 SAMPLES/CYCLE 
A - 1 (KNOWN) 

0.7 l I I I I I — 1 I I I I I I I 1_ I 1 I - 

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1718 

NUMBER OF SAMPLES 


Fig. 2. Least squares algorithm: noise free case and 
baseband sampling, n = 5 










